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Abstract. We present global existence results for solutions of reaction-diffusion systems on 
evolving domains. Global existence results for a class of reaction-diffusion systems on fixed 
domains are extended to the same systems posed on spatially linear isotropically evolving 
domains. The results hold without any assumptions on the sign of the growth rate. The 
analysis is valid for many systems that commonly arise in the theory of pattern formation. 
We present numerical results illustrating our theoretical findings. 



1. Introduction 



Since their seminal introduction by Turing [1952], reaction-diffusion systems (RDS's) have con- 



stituted a standard framework for the mathematical modelling of pattern formation in chemistry 
and biology. Recent advances in mathematical modelling and developmental biology identify 
the important role of domain evolution as central in the formation of patterns, both empiri- 
cally Kondo and Asai 1995 and computationally Comanici and Golubitsky 2008 Crampin 



et al. 1999 Madzvamuse and Maini 



Madzvamuse! 2006 and Barrass et al 



2007 . In this respect, many numerical studies, such as 



2006 , of RDS's on evolving domains are available. Yet, 



fundamental mathematical questions such as existence and regularity of solutions of RDS's on 



evolving domains remains an important open question Kelkel and Surulescu 2009 



Numerous studies on the stability of solutions of RDS's on fixed domains are available, for 
example, iHollis et all 119871; iRothel 11984]; |Wei and Winterl 12008] , but very little literature 



regarding the stability of solutions of RDS's on evolving domains. Madzvamuse et al. 2010 
provides a linear stability analysis of RDS's on continuously evolving domains, and Labadie| 
2008 examines the stability of solutions of RDS's on monotonically growing surfaces 



Oi 



discussion here differs from all these studies in that we focus on planar evolving domains and 
we show existence, uniqueness and stability, for an entire class of RDS's on evolving domains 
independently of the rate of evolution. In this article we prove the stability of solutions of RDS's 
on a particular, but fundamentally important, class of time-evolving domains: that of bounded 
spatially linear isotropically evolving domains. 

We show that if a RDS fulfils a restricted version of certain stability conditions, introduced 
by Morgan 1989 for fixed domains, then the RDS fulfils the same stability conditions on any 



bounded spatially linear isotropic evolution of the domain. We thus prove that, under certain 
conditions, the existence and uniqueness for a RDS on a fixed domain implies the existence and 
uniqueness for the corresponding RDS on an evolving domain. This is, to our best knowledge, 
the first result that holds independently of the growth rate and is thus valid on growing or 
contracting domains as well as domains that exhibit periods of growth and periods of contrac- 
tion. Our analysis rigourously justifies computations for this type of domain evolution, and we 
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illustrate our results with benchmark computations using a moving finite element (Lagrangian) 
approach. 

The outline of our discussion goes as follows: In ^2]we state our model problem together with the 
form of domain evolution that we consider, and present a transformation of our model system 
to the Lagrangian framework that helps in proving global existence of solutions. In plwe review 
the existence results for RDS's on fixed domains which will form the basis of our analysis. In £J4] 
we state and prove the central results of this work; in particular, we extend the existence results 
cited in |3]to problems posed on evolving domains. In S|5]we illustrate some specific applications 
of our results, in particular those of significance in the field of biological pattern formation. We 
focus on growth functions commonly encountered in the field of developmental biology for which 
our analysis is valid and show the applicability of our analysis to some of the important reaction 
kinetics encountered in the theory of biological pattern formation. In ^6] we present numerical 
results for a RDS posed on a periodically evolving domain. We present a moving finite clement 
scheme and a fixed domain finite element scheme to approximate the solution of a RDS posed 
on the evolving and the Lagrangian frame respectively. In Sj7] we summarise our findings and 
indicate future research directions. 



2. RDS'S ON CONTINUOUSLY EVOLVING DOMAINS 

Let u(x,t) be a (m x 1) vector of concentrations of chemical species, with x £ K ra , the time- 
dependent spatial variable and t £ [0,T], T > 0, the time variable. The model problem we 



wish to consider is a semilinear RDS posed on a continuously evolving domain (see Madzvamuse 
1 2000 1 for details of the derivation), given by, 

(2.1) d t u (x, t) + [V • (a : u)} (x, t) = f(u (x, t)) + DAu (x, t) for x £ Q t and t £ (0, T], 

where fl t C K™ ( n < oo) is a C 2+7 (S1), simply connected, bounded and continuously deforming 
domain with respect to t. The function / is a (m x 1) vector of nonlinear coupling terms that 
is locally Lipschitz, D is a diagonal matrix with strictly positive entries on the diagonal, a is 
a flow velocity generated by the evolution of the domain and (a : u) := (au\, . . . ,au m ) T by 
definition. 

2.1. Assumption (Flow velocity) We assume that the flow velocity a,i(x,t) is identical to the 
domain velocity, i.e., 

a, = d t x % i = 1, ...,n, 
as is standard in the derivation of RDS's on evolving domains on application of Reynold's Trans- 



port Theorem Acheson 1990 



To simplify the exposition we take boundary conditions to be of homogenous Neumann type. 



Morgan 1989 considers a much wider class of boundary conditions and our analysis may be 
extended to this more general setting. We are primarily interested in patterns that arise as 
a result of self-organisation prompting the consideration of homogenous Neumann boundary 
conditions. We take the initial condition for each Ui to be bounded and nonnegative. Our 



model problem thus takes the form (2.1 1, equipped with the following boundary and initial 
conditions: 

f>-vu](z,i) = o, X £dn t ,t>o, 

1 u(x,0) = u (x), x£fl . 



GLOBAL EXISTENCE FOR SEMILINEAR REACTION-DIFFUSION SYSTEMS ON EVOLVING DOMAINS 3 



2.2. Lagrangian transformation. For our analysis it is more convenient to work with problems 
denned on a time-independent domain. Therefore, we introduce a transformation that maps our 



model problem (2.1) from a time-dependent domain to a fixed domain (see Baines 1994 for 



a more detailed discussion of this approach). In order to do this, without too many techni- 
cal complications, we will restrict our attention to special evolutions of the domain, described 
next. 

2.3. Assumption (Isotropic domain evolution) We assume the domain H, t to evolve by obeying 
a bounded spatially linear isotropic domain deformation, i.e., 

(2.3) x = p(t)£ for all (£, t) £ fl Q x [0, T] and all x £ Sl t , 

with p £ C 2 ([0, T]; 0, oo) and where £ represents the spatial coordinates of the initial domain. 
This assumption and Assumption |2 . 1 1 imply that 



(2.4) 
wher 



a(x,t) = />(*)£, 



e p 



._ dp 
dt ' 



Hence, we obtain the following transformed problem on the initial domain (see Madzvamuse and 
Maini| |2007) for details), with 

for ie [0,T] and£e fi , 



(2.5) 
we have 

(2.6) 



u(£,t) = u(p{t)£,t) 



dtii + n-u 



[p-vu]ie,t) = o 

a(€.o) = fio(€) 

< «o(0 < oo, 



p 2 



onO x (0,T], 
on<9O ,t>0, 



where Slo is the initial spatial domain, n is the spatial dimension and the Laplacian is now 
taken with respect to £. The local and global existence results that we utilise from the existing 
literature require the coefficients on our diffusion term to be independent of time; to this end we 
introduce the following proposition: 

2.4. Proposition (Time rescaling Labadie 2008 ) Let u be a solution of (2.1|, rescaling time 
via the change of variables 



(2.7) 



»(*) := 



d ■;■ 



/o P( r ) 2 ' 
and denoting S := s(T). We have, u (p(t)£,t) = u(£, s), where u satisfies 

{d s u + np'pu = p 2 f(u) + DAu 
p.Vtt](&*) = 
«(€,o) = t«o(€) 



onfi x (0,S], 
££dn,s>0, 



Furthermore, if f(u) is locally Lipschitz in u then f(u(£,, s), s) — p 2 (s)f(u(£, s))—np(s) p(s)u(£, s), 
is locally Lipschitz in u. 

Proof We note that with domain evolution of the form considered in this study, there exist 



C\,C<i such that < p < C\ < oo and HpIIl^^t] < C2 < oo. Applying the rescaling (2.7|, we 
see that for any function g £ C 1 [0, T] 

(2.9) d s g (t) = d s t(s)d t g(t) = p 2 (t)d t g{t). 
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Defining m(£. s) :— u(p(t)£,t) and multiplying problem (2.6) by p , we obtain problem (2.8) 



Clearly since f(u) is locally Lipschitz in u, Assumption 2.3 implies f(u, s) is locally Lipschitz 
in u and globally Lipschitz in s. □ 



3. Basic theoretical setting for RDS's on fixed domains 



We now summarise the existing results for RDS's on fixed domains which will form the basis of 
our analysis on evolving domains. The following result is a straightforward generalisation of 
Hollis et al.| [l987| Prop. 1]. 



3.1. Theorem (Local existence) Let Assumption 2.3 hold. Problem (2.8) admits a unique local 
solution. Furthermore, defining the unique maximal solution of (2.8) by 



(3.1) u:O x[0,T max )^M- 
there exists a function TV € C([0,T max );IR m ) such that, 

(3.2) \\ui(;s)\\ Loo (n o) <Ni(s) forie [l,...,m] and s G [0,T max ). 
Finally, if T max < oo, 



(3.3) 



lim 



£IM 



•iSJIUoo(n ) 



i=l 



3.2. Global existence on fixed domains. If p(t) = 1 for all t £ [0, T], Problem (2.8) becomes, 
find u-M Q x (0, T] -> W n such that, 



(3.4) 



d t u = /(«) + DAm, on Ct X (0, T], 
[*7-Vi2](£,*)=0, £ed£W>0, 

u(C,0) = iio(O. ^e^o, 



where we have used the fact that s(t) — t (cf. ( |2.7| ). 

3.3. Definition (Invariant region) S C M m is called an invariant region for the solution of the 
reaction-diffusion system (3.4) if for any solution u, 



(3.5) «(£,0) e S =» u(£,t) e S for all te (0,T]. 
3.4. Assumption (Positive solutions) We assume hereon that 

(3.6) /»(«)k=o > for all te [0,T], 

and for / ^ C^M^jM 1 "), the strict inequality 



(3.7) 



/,-Mk=o > for all te [0,T]. 



Assumption 3.4 together with the positivity of our initial data, implies WV which we refer to as 



the positive quadrant, is an invariant region for the solutions of problem (3.4) (see Smoller 1994 
Th.14.7,14.11 pp.200-203]). 



3.5. Remark (General invariant regions) Assumption 3.4 may be relaxed. The proof of our 



existence results only requires bounded initial data and the existence of an invariant region. 
Furthermore, consideration of the positive quadrant alone is sufficient for our studies. 
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3.6. Lyapunov stability conditions. We now introduce a Lyapunov function for the dynamical 
system defined by ( 3.4 ) when the initial condition Uq varies which is used to prove global existence 



and a restricted version of the conditions it is required to fulfil Morgan 1989 



Suppose / is as defined in problem (3.4) and that there exists a function H E C 2 (R + ;R) and 
hi G C 2 (R+; R) for each i = 1, . . . , m, such that 



(3.8) 



H(z) = S~] hi(Zi) for all z <G 



(3.9) hi(Zi),h%(Zi) > for all z e R™ 

(3.10) H(z) -^oo <^=^ z^oo for all z£R™. 

Suppose there exists A = (oy) e (R) mxm satisfying ay > 0, an > with 1 < i,j < m such that 
for some r, k\, k% € R+ independent of j, we have 

3 

(3.11) ^Vy/i^)/^) < ki(H(z)) r + k 2 for all z e R™, j < m. 

Suppose there exist q, k^,ki G R+ such that for 1 < i < m, we have 

(3.12) / l , i (^)/ i (z)<fe 3 (F(^))«+fc4, for all zeE™ 
Suppose there exist k§ , fcg > such that 

(3.13) Vtf(*) • /(*) < k 5 H (z) + k 6 for all z e R™. 

3.7. Theorem (A priori estimates (Morgan! ' 1989 1) Let conditions (3.8), (3.9) and (3.13) hold 
and let it be a solution of problem (3.4). The following a priori estimates hold, 

rt 

<g(t) forO<T<i<T max , 



(3.14) 
(3.15) 



Loc (Ho) 



J H(u(; S ))ds 

j j H{u(£,s)) 2 d£ds<g{t) iorO<t<T n - 



where g,g € C[0, oo). 

3.8. Theorem (Global existence on fixed domains [Morg an 1989 ) If conditions (3.8) — (3.12) 
hold, with r from condition ( 3.11 1 satisfying r < (1 + a), a £ R + , u is a solution of problem (3.4 ) 
and if there exists g 6 C[0, oo) such that 
ft 



(3.16) 



H(U(; S )) 



ioo(0 ) 



< g{t) for < t < t < T a 



then T max = oo. Alternatively, if conditions ( |3.8[ ) — (3.12) hold, with r from condition (3.11) 
satisfying r < (1 + ^zr^), b > and where n represents the spatial dimension, u solves a problem 
of the form (|3.4| and if there exists g € C[0, oo) such that 



(3.17) 



then Trr 



/ ( H{u(i,s)) d£ds<g(t) for < t < T max , 

JO JSln 



Specifically if r from condition (3.11 1 satisfies r < 2 or if Slo C R with r < | and the remaining 
conditions (3.8) — (3.13) are satisfied then T max = oo. 
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4. Global existence on evolving domains 



In this section we show that, if the stability conditions in [3.6 are valid for Problem (3.4) then 



they remain valid under any evolution of the domain fulfilling Assumption |2.3| given a suitable 
assumption on the structure of H. We also extend the previous a priori estimates and existence 
results of Morganj [1989 to problems with time dependent /. 



4.1. Assumption (Polynomial Lyapunov function) We assume the Lyapunov function introduced 
in §3.6| is of the following form 



(4.1) 



H{z) = Y j zf, Pi >liovi = l, 



. ,m. 



s=l 



4.2. Remark (Polynomial growth restriction) Assumption 4.1 is somewhat natural. Condition 



(3.12 1 is essentially a polynomial type growth restriction on the zero order terms Morgan 1989 



Assumption 4.1 can be viewed as the explicit analogue of the polynomial growth restriction on 



the zero order terms implicit in (3.12) 



4.3. Lemma (Equivalence of Lyapunov functions) Suppose Assumptions 2.3 3.4 and 4.1 hold. 

Let t he Lyapun ov sta bility conditions in §3.6| be satisfied by H and /. Then the con ditions in 

3.6 with r (cf. (3.11)) replaced by f := max(l,r), are satisfied by H and / (cf. (4.2)) in place 



of/. 



Proof We denote the zero order term in Problem (2.8) by 

(4.2) / (u{l a), a) := p 2 {s)f (u{l a)) - np{s)p{s)u{l s). 

The positive quadrant remains an invariant region for the solutions of our evolving domain 
problem since 



(4.3) 



fi [u(C s),s) |„ i=0 = p 2 {s)fi (u(i,a)) | Ui=0 , 



problem (3.4). Denote C\ :- 



\\P\ 



thus Assumption 3.4 implies M™ is an invariant region for the solutions of problem (2.8). Let 
ki,i = 1, ...,6, q,r and A be as defined in [3.6 for which conditions (3.8) — (3.13) hold for 

IHIz,oo[0,T] and Ca 



£oo[0,T] 



; these arc well defined real 



numbers thanks to to Assumption 2.3 We now show that conditions (3.8) — (3.13) hold with the 



same H, f replaced by / and r from (3.11 ) replaced by f, where f — max(l, r). 



Clearly conditions (3.8) — (3.10) are still satisfied as they depend only on H which is unchanged. 
Condition ( |3.11[ ) holds since 

3 i i 

(4.4) J2 a v h 'ih = Y, a iMp 2 fi - n PP^i) < (ki(H) r + k 2 )Cl + nCiC 2 £ cHjhfa, 

i=l i=l 

by the stability of the fixed domain problem. Assumption |4.1| gives, 
i j 

^a^h'Ji < (h(HY + k 2 )C 2 + ndC^aimhi 



r-i 



(4.5) 
where 



8=1 

^2 i i„ tj ^ fi. /^i2 



< (fci(# ) r + k 2 )C( + k 7 H < (kiCf + k 7 )(H) r + k 8 , 



max(l,r) and k s := < 



k 2 C\ + hCl ifr<l, 
k 2 Cl ifr = l, 

k 2 Cl + k 7l ifr>l. 
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Condition (3.12| holds since 

Kfi^Clh'Ji + nCxCiPihi 

(4.6) < hCf(H) q + k 4 Cl + ndC 2 max( Pl )H 

<fcioff 9 ~ + fcii, 



where 



k±Cf + k 3 Cl if g< 1, 

q = max(l, q) and fen := ( k^Cf, if q = 1, 

fc-iCf + fc 9 , if g > 1. 



Condition (3.131 holds since 

VE/ < Y, C l h 'ifi + nC 1 C 2 p i h i < (k 5 H + k 6 )C* + k 12 H. 



(4.7) 



i=l 



Thus the positive quadrant remains an invariant region for the solutions of problem ( 2 
the Lyapunov stability conditions in §3.6|are satisfied, completing the proof. 



and 



□ 



4.4. Remark (Applicability of Morgan |1989 to systems with time dependent zero order terms) 
Suppose the reaction function / (u(£,t),t) is locally Lipschitz with respect to u and t, and 
suppose that the Lyapunov function H depends only on u. Then, Theorems 3.7 and 3.8 remain 



applicable |Morgan 1989 (5.5)], jMorgan and Holhs 1995] Th. 1.1] and Bendahmane and Saad 



2010 Th. 4]. Thus, the Lipschitz result of Proposition |2.4| the structural Assumption 4.1 and 



the equivalence of Lyapunov functions proved in Lemma |4.3| imply that Theorems |3.7| and |3.8 



are applicable for solutions of (2.8 1. 



For completeness, we include a proof of Theorems 3.7 and 3.8 for solutions of Problem (2.8 1, in 



Appendices [A] and IB] respectively. To remain concise we prove a sufficient existence result for 
the examples presented in Efs] and briefly sketch the full proof of Theorems |3.7| and |3. 



The results of Morgan and Hollis 1995 apply for systems with time dependent diffusion. This 
may allow treatment of more general domain evolution where the rescaling carried out in 32] 
yields a system with time dependent diffusion. We leave this generalisation for future studies. 



4.5. Theorem (Global existence of solutions on evolving d oma ins) Let Assumptions 2.3 3.4 and 

,r<§ (cf. 



4.1 



hold and suppose H, f and r satisfy the conditions in 1 3.6 with r < 2 or if Q C 



(3.11)). Then, Problem (|2.l[) admits a global classical solution. 



Proof Application of the results in § 2. 2 1 allows us to show existence for the transformed Problem 



classical solution. From Lemma 



and r < 2 or if f2 C 



4.1 



(2.8) defined on a fixed domain. Theorem 3.1 gives the existence of a unique non-continuable 



r < 



Theorem 



3.7 



t he st ability conditions in [3.6 hold with / (cf. (4.2)), H 



gives an a priori estimate for H. Theorem 3.8 implies 



?max = co (cf. (3.1)) completing the proof. 



□ 



5. Applications 



In this section we illustrate some applications of Theorem |4.5[ We present different forms of 
admissible domain evolution that fufil Assumption |2.3| We show that Assumption |4.1| is applica- 
ble to some commonly encountered models in chemistry and biology. We concentrate on RDS's 
which admit Turing instabilities, as the main focus of our research is biological pattern formation. 
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We identify and describe Lyapunov functions and constants that imply global existence for the 



fixed domain problem and thus for the evolving domain problem by Theorem 4.5 



5.1. Admissible domain evolution. We now provide some commonly encountered examples 
of domain evolution in developmental biology for which Assumption |2.3| holds: 



• Logistic evolution on any finite positive time interval 



(5.1) 



P(t) = 



l + i( e r,t_l) 



, t£[0,T], 



where r g > is the growth rate and K > 1 is the carrying capacity (limiting size of the 
evolving domain). 

• Exponential evolution on any finite positive time interval 

(5.2) P(t)=e r »*, te[0,T}. 

• Linear evolution on any finite positive time interval 

(5.3) p(t) = l + r g t, te[0,T], 
where r g > —yp. 



5.2. Admissible kinetics. We now present some of the commonly encountered reaction kinetics 



of problem ( 2.1 1 for which the analysis of Morgan 



1989 implies global existence of solutions on 



fixed domains. We first consider the problem (|2.l|) with this general reaction term 

(5.4) 



Iii u ) = ^2djUj + (-l) l g(u) + h, 

3=1 



with the following restrictions 

(5.5) 
(5.6) 
(5.7) 



Cij > 0, for i 7^ j. 



(5.8) 



(5.9) 



bi>0, i = l,.. 

/(■u)U i= o = 0, i = l,.. 



, m. 
, m. 



for all u <G 



•■+• 



The motivation of this type of kinetics, as discussed by Murray 2003 , is their role in the theory 
of biological oscillators due to a feedback mechanism. 



5.3. Proposition (Lyapunov function) We show that problem (2.1 ) equipped with kinetics (5.4) 
is well posed, with Lyapunov function H{z) := Y^hLi z i- 



Proof Recalling that the initial data is bounded and nonnegative (2.2), we show that Assump- 
tion 3.4 is fulfilled, which implies E+ is an invariant region for the solutions. Indeed, from (5.7) 



we have 
(5.10) 



fi(u)\ Ui=0 = y^ cjjUj 

3 = 1 



i—1 rn 
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Conditions (5.5) and (5.6) imply 

(5.11) fi(u)\ Ui = >h>Q for all ue 



with r — 1. Clearly conditions (3.8) — (3.10) hold. Condition (3.11) holds with 



Thus Assumption 3.4 is fulfilled due to (5.9). Now we show conditions (3.8) — (3.13) are fulfilled 



1 if j = i, 

1 if j = 1 and i is even, 

otherwise, 



fci = maxij(cy), &2 — 5Zi=i fr« anc ^ r ~ 1- Condition (3.12) holds since 

ni 
(5.12) /i'/j = fi< max(ci,) Y^ m + max(6, 4 ) + g(u). 



1=1 



Using (5.8) we have 
(5.13) 



M/i < max(c„) ^ Uj + max(6i) + (} Uif 

i,j * — * i f — ' 



/ = J i=l 



< max(c y ;)-H"(m) + i/(u) p + max(6,) < fciJ(u) 



'■./ 



where k,q,b G 
holds since 

(5.14) 



represent constants that depend on the value of p in (5.8). Condition (3.13) 



VH-f = J2f*< max(cy) Vui+Vl»i< fciJ(w) + 6, 



i=l 



i=l 



where fc = max,;j(cjj) and b — Y]j—i bj. A straightforward application of Theorem 
completes the proof. 



4.5 



□ 



thus 



5.4. Examples. The generic problem for which we showed global existence of solutions actually 
encompasses some of the more widely studied models in the theory of pattern formation such 
as the Gray-Scott model and the Brussellator. Below we present two examples of two species 
reaction terms for which our analysis implies global existence of solutions, the first of which is a 
restriction of the reaction term above. 

• Activator- depleted substrate model: We consider the activator-depleted substrate model 



Gierer and Meinhardt, 1972 Lefever and Prigogine 1968 Schnakenberg 1979 also 



(5.15) 



known as the Brusselator model: 

\fi{ui,u 2 ) =i(a-ului), 
\/2(-"i,u 2 ) = 7 (& - «2 + wfwi) , 

where < a, b, 7 < 00. The assumption of nonnegative initial data implies that the 
positive quadrant is invariant for our problem due to the fact that a, b > 0. If we take 
H(u) = U1+U2, then conditions (3.8) — (3.13) are fulfilled with r — which implies global 



existence of solutions on evolving domains via Theorem |4.5| The remaining constants 



for which conditions (3.8) — (3.13) hold are given in Table [T 



• Thomas reaction kinetics: The following model, proposed and studied experimentally by 



Thomas 1975 , is based on a specific reaction involving the substrates oxygen and uric 



10 
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(5.16) 



(5.17) 



acid which react in the presence of the enzyme uricase: 

\fi(ui,u 2 ) =7(a-ui- g(u 1 ,u 2 )) , 
\f2(ui,u 2 ) = j(b-au 2 -g(u 1 ,u 2 )), 



where 



g(ui,u 2 ) 



KU\U 2 



1 



and < 7, a, a, b, K, /3 < oo. Once again the assumption of nonnegative initial data 
implies that the positive quadrant is invariant for our problem due to the fact that 
a, b > 0. If we again take H(u) — u% +u 2 , then conditions (3.8| — (3.13) are fulfilled with 



r = which implies global existence of solutions on evolving domains via Theorem |4.5| 
The remaining constants for which conditions (3.8) — (|3.13l hold are given in Table nl 



5.5. Remark (Invariant rectangles for the Thomas model) It can be shown, utilising the tech- 
niques of Smoller |1994| , that there exist bounded invariant rectangles for the solutions of the 
Thomas model defined in (5.16). This implies global existence of solutions via Theorem 3.1 



However the authors can show that it is possible to construct growth functions that fulfil As- 
sumption 2.3 for which the bounded invariant rectangle can be made arbitrarily large. This 



necessitates the Lyapunov function approach to show existence and uniqueness of solutions. 



Parameters 
an 
ai2 

021 

022 
fcl 

A'2 

k 3 

<1 

h 



Activator-depleted substrate model 

1 



1 

1 



j(a + 
1 

2 

3 



•6) 



7(max(a, b)) 



7(0 + b) 



defined in (5.15) and (5.16) respectively. 



Table 1. Terms from {3 for which conditions (3.8) — (3.13) hold for the kinetics 



Thomas model 

1 





1 


7(0 + b) 




7(max(a, b)) 


j(a + b) 



5.6. Remark (Further applications) The analysis can be applied to a large number of problems 
unrelated to the theory of pattern formation. Garvie and Trencheal 2009 provide an example 



applicable to ecology and the aforementioned paper of Morgan 1989) contains further examples 



as well as the numerous citations of said paper that use the approach on various problems. 



6. Numerical experiments 



In this section we present numerical results on two-dimensional evolving domains to back-up 
the theoretical results of the previous sections. The numerical simulation of RDS's on growing 



domains is an extensive research area. Crampin et al. 2002 and Madzvamuse and Maini 2007 



study mode doubling and tripling behaviour of the activator-depleted substrate model on one- and 
two-dimensional growing domains. The theoretical results derived above apply for any evolution 
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that fulfils Assumption 2.3 We present numerical results on a periodically evolving domain that 
exhibit spot splitting as well as spot annihilation and merging, which to the authors knowledge 
is as yet an unstudied area. 

6.f. Domain evolution. We consider periodic domain evolution defined by 

(6.1) p(i) = l + 9sin(^J , t € [0,1000 = T], 

with the initial domain defined as Slo = [— 0.25,0. 25] 2 which grows to fJsoo = [— 2.5,2. 5] 2 before 
contracting back to the original domain. 

6.2. Continuous problems. We present results for the aforementioned activator-depleted model 
considered on a periodically evolving domain and its equivalent transformed system on a fixed 



domain as in { 2.2 For example, on a periodically evolving domain the problem is stated as 

follows: 

(6.2) 

f 

d t ui(x,t) + [V • (aui)](x,t) — Aui(x,t) = 0.1 — [ufwiKa^t), for x G Q t 

d t u 2 (x,t) + [V- (au2)](x,t)-0.01Au2(a:,t) = 0-9 - u 2 {x,t) + [v%ui](x,t), and t G (0,T], 
[v- Vu](x,t) = 0, xedn t ,t>0, 

where x = p(t)£ and a,i(x,i) — p(t)£, i = 1,2. Equivalently, the following transformed equations 
are obtained on a fixed domain, 

d t ui + 2fui - jjAm = 0.1 - u\u Xl 
Hi.3) {d t u 2 + 2^u 2 -^Au 2 =0.9-u 2 + ulu 1 , onfi x(0,T], 

i?-Vu](£,t)=0, £edn a ,t>0. 

In both cases we take identical initial conditions as small perturbations around the homogenous 
steady state of (1.0, 0.9) obtained in the absence of domain evolution. 

6.3. Numerical schemes. We employ a Galerkin finite element method for the spatial approx- 
imation and an implicit-explicit modified backward Euler scheme for the time integration. Dis- 
cretising in time we divide the time interval [0, T] into a partition of N uniform subintervals, 
= to < . . . < ijy = T and denote by r := t n — t n —i the time step. For the spatial discretisation 
we introduce a regular triangulation T° of f^o with K G T° an open simplex. We define the 
following shorthand for a function of time, f(t n ) ='■ /". 

We define the finite clement space on the initial domain V° C H 1 (Q,q) as, 

(6.4) V° := {$ G H\too) : ®\ K € P 1 for all K G 7^}, 

where P 1 denotes the space of polynomials no higher than degree 1. For the numerical simulation 



of equation (6.2) we require finite element spaces defined on the evolving domain. We construct 
the finite element spaces V™ according to the following relation between the basis functions of 
V" and V°. 

(6.5) * n = *(/t> n €, *„) = *«) n=l,...,N. 
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Thus the family of finite element spaces on the evolving domain V" C H 1 (fi t?i ) n 
may be defined as, 



,N 



(6.6) 



:= {*" e H^fltJ : $ n \ K e P 1 for all K e T™}, 



where we have used the fact that the domain evolution is linear with respect to space. We 
approximate the initial conditions in both schemes by 

(6.7) U° = lu {£) for all £ € Oq, 

where I is the standard Lagrange interpolant. The finite element scheme to approximate the 



solution to equation (6.2 ) aims to find U{ 1 , U% € V", n = 1, . . . , N such that 




i(kT,*") + (Wf, W") = (0.1 - (u^ 1 ) U{\W l ) + i^™- 1 ,^"- 1 ), 

±{U$, *"> + 0.01<W.™ V*") = (0.9 - [/« + U^UfUg, #") + ^(U^' 1 ,^ 71 ' 1 ) 



, n 



1, . . . , N. Similarly the finite element scheme to approximate the solution 
to equation (6.3 ) aims to find W™, Wg £ V°, n = 1, . . . , N such that 



wj* -1 ,*) 



1(W? 
for all $ € V°. 



^(vwf, v$) + ^(wr,*) = (o.i - ( ir : r 



W?,*>, 



^(VW 2 ",V$) 



2^; 



(Wa™, $> = (0.9 - W£ + W^-'W^W?, $) 



We solved the models in C utilising the FEM library ALBERTA by Schmidt and Siebert 2005 



We used the conjugate gradient solver to compute our discrete solutions. We took an initial 
triangulation T° with 8321 nodes, a uniform mesh diameter of 2~ 6 and a fixed timestep of 10~ 3 . 
PARAVIEW was used to display our results. 



6.4. Results. Figures [T] a nd [2| show snapshots of the activator profile corresponding to the 
activator-depleted system (5.15). The inhibitor profiles have been omitted as they are 180° 



out of phase to the activator profiles. We have verified numerically that there is very little dif- 



ference between the discrete solution corresponding to system (6.8| mapped to the fixed domain 



and the discrete solution corresponding to system (6.9) defined on a fixed domain, as is expected 
from the results in 92.21 



The figures illustrate the mode doubling phenomena that occurs as the domain grows as well 
as the spot annihilation and spot merging phenomena that occurs as the domain contracts. In 
Figure [3] we present in more detail the novel spot merging phenomena observed on the contracting 
domain. It is still unclear whether the spot merging phenomenon is in fact a special case of the 
spot annihilation phenomenon, that occurs when the modes are of sufficient proximity to influence 
each other. 

We note that the mode transition sequence, i.e., the number of spots, is different when the 
domain grows to when it contracts. The difference in the mechanism of mode transitions on 
growing and contracting domains is an area in which very little work has been done and these 
initial numerical results indicate the need for further exploration of this area. 
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7. Conclusion 

Many problems in biology and biomedicine involve growth. In developmental biology recent ad- 
vances in experimental data collection allow experimentalists to capture the emergence of pattern 
structure formation during growth development of the organism or species. Such experiments 
include the formation of spot patterns on the surface of the eel, patterns emerging on the surface 
of the Japanese flounder and butterfly wing patterns forming during the growth development of 
the imaginal wing disc. In all these examples, patterns form during growth development. 



Since the seminal paper by Turing [1952| which considered linear models that could give rise 
to spatiotemporal solutions on fixed domains due the process of diffusion-driven instability, a 
lot of theoretical results on global existence of such solutions have been derived and proved for 
highly nonlinear mathematical models Rothe 1984 ; Smoller 1994 . Only recently, mathematical 



models on growing domains have been derived from first principles in order to incorporate the 



effects of domain evolution into the models Crampin et al. 1999 ; Madzvamuse 2000 . In all these 



studies, very little analysis has been done up to now to extend the theoretical global existence 
results to models defined on evolving domains. 

Under suitable assumptions, we have extended existence results from problems posed on fixed 
domains to problems posed on an evolving domain. We have illustrated the applicability of 
the existence results of Morgan 1989 to problems on evolving domains. We have shown that 



global existence of solutions to many commonly encountered RDS's on fixed domains implies 
global existence of solutions to the same RDS's on a class of evolving domains. The results are 
significant in the theory of pattern formation especially in fields such as developmental biology 
where problems posed on evolving domains are commonly encountered. Our results hold with 
no assumptions on the sign of the growth rate, which may prove useful in other fields where 
monotonic domain growth is not valid from a modelling perspective. The applicability of our 
results is demonstrated by considering different forms of domain evolution (linear, logistic and 
exponential) . 

In order to validate our theoretical findings, we presented results on a periodically evolving do- 
main. Our results illustrate the well-known period-doubling phenomenon during domain growth 
but more interesting and surprising is the development of spot annihilation and spot merging 
phenomena during contraction. This raises new questions about bifurcation analysis on growing 
and contracting domains. 

One of our primary goals is the numerical analysis of finite element approximations of RDS's on 
evolving domains. The classical existence results obtained will be an important tool in future 
work. Numerical experiments have been carried out and they illustrate the need for further 
numerical analysis especially in the case of contracting domains. Extension of our work onto 
domains with more complex evolution is another area for future research. 
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Figure 1. Snapshots of the discrete solution U% corresponding to system (6. 



at times 0, 50, 160, 220, 380, 500, 700, 740, 820, 900, 980 and 1000 reading from 
left to right and then top to bottom. For parameter and numerical values see 
§6.3| The solution exhibits a mode doubling sequence of 1,2,4, 8 and finally 10 
as the domain grows. As the domain contracts the spots are annihilated in a 
sequence of 8, 6, 4, 2 and the final transition to a single spot occurs via merging, 
with the final domain exhibiting no patterns. 
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Figure 2. Snapshots of the discrete solution W2 corresponding to system (6.9 1 
at times 0, 50, 160, 220, 380, 500, 700, 740, 820, 900, 980 and 1000 reading from 
left to right and then top to bottom. For parameter and numerical values see 



§6.3| The mode transition follows exactly that of Figure [TJ corroborating the 
results in 92.2 
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Figure 3. Snapshots of the discrete solution corresponding to system (6.8) 



(top) and system (|6T9j) (bottom) at times 930, 940, 950, 955, 965, 975. The spot 
merging phenomena observed in the transition from two spots to one spot is 
displayed. 
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Appendix A. A priori estimates for systems with time dependent / 



In this section we prove Theorem 3.7 remains applicable for solutions of Problem (12.81). For 



conciseness we focus on assertion (3.141. We adapt the proof of Morgan 1989 Th. 3.2] to our 
purposes. In Appendix [B] we use the a priori estimates obtained in this section to prove the 
global existence of classical solutions to Problem (2.8 1. We state the main result of this section 
in the following Lemma. 



A.l. Lemma (An a priori estimate for solutions of Problem (2.8)) Suppose Assumptions 2.3 



3.4 and 4.1 hold and let (3.8) — (3.10) and (3.13) hold. Then, assertion (3.14) is valid with u a 



solution of Problem (2.8). 



Proof Let / be as defined in (4.2). From the proof of Lemma 
remainder of the proof into steps. 

Step 1: We first show the following inequality for H: 



4.1 



(4.7) holds. We split the 



(A.1) 



(A.2) 



J T A — 1 



+ H (u{t r)J + k 6 Cf(t - t) for £ e n and t < T max 
From J2J8) we have for (£,*) e Q X (0,T max ), 

VH («(£*)) ■d t u(g,t) = VH(u(£, t))-(p Au(£,t) + /(«(£, t),t) 



Using ( |3.8[ ), ( |4.7[ ) and (A.2) we obtain the following generalisation of Morgan 1989 
(3.1)], 



(A.3) VH («(£, t)\ ■ d t u(£, t)<J2 AA^ (u;(f, t)\ + (fcCf + k l2 ) H (u{£, t)\ + k^G\ 



where we have used the convexity of H (3.9). Integrating (A.3 1 in time gives (A.l). 



Step 2: We use (A.l) to construct an appropriate barrier function with a view to applying the 
maximum principle, corresponding to Morgan 1989 (3.2) — (3.4)]. Introducing an arbi- 
trary T* < T max , we define 

/■t rn, 



(A.4) 



D, 



Y J^hi (ui(£,r))dr for £ e fl and < r < t < T* , 



(A.5) 



where D* := max(_D^). Observe that from (3. 

i 
J H(u(; S ))ds 



and (A.4 1, we have for < r < t < T* 



£00(^0) 



r _ _ 
/ ^fti(ui(-,s))ds 



< max 



D* 



l|w(-,s)ds|| Loo(no) . 



From (A.l) and (|A.4[) we obtain 



(A.6) 



ftu;(€, t) < D*Aw(£, t)+L + Mwfa t) + k & Cl{t - r), ?eJ](,,0<r<t<r 



[i/-VH(€.*)=0, 
w(€,t)=0, 



^ai]„,te(r,r] 



IS 



(A.7) 
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where 



L:=\\H(u(;t)) lU^no) and M := 



(fc 12 + Clh) D* 
min (Di) 



(A. 



(A.9) 



Note we have used ( 2.8 ) to obtain the boundary conditions. For the purposes of applying 



the maximum principle we define a barrier function 

w(£, t) ■- w(£, t) - L + K e C i T for £* e Q and < t < t < T* 



From (A.6) and (A.4) we have 



d t w(& t) < D*Aw(£, t) + Mw{i, t) £en Q ,0<T<t<T* 

[i/-v«>](|,t) = o, ^3n ,te(T,r] 

w(Z,t)<0, £en . 



Step 3: We use the maximum principle to complete the proof. Applying the strong maximum 



principle for parabolic problems Sperb 1981 Th. 2.9, Rem. (a) pg. 21] to (A.9) and 
noting the positivity of w, we have 



(A.10) 



(A.11) 



(A.12) 



L + K e CfT* 
M 



< w{£, t) < for all £ € n Q and for t e (r, T*]. 



From (A.8) and (A. 10) we have 



< w(£,t) < L + KeC i T f or a n ^ G q q and for t e ( T ,T*]. 



We conclude from (A. 11) that 



fit" (~l2rp* 

\\w(;t)U x( n o) < L + K ^ 11 fate(r,r]. 



Since T* was arbitrary, combining (A. 5) and (A.12) completes the proof of the Lemma. 

□ 



For completeness, we sketch the proof of assertion (3.15) with u a solution of Problem (2.8). In 
(A.3) we denote K 7 := k 5 Cf + k i2 and K s := k 6 C(, where K 7 ,K S correspond to the terms on 



the right hand side of |Morgan 1989| (3.1)]. Assertion (3.15) follows from the proofs of Morgan 
T989] Th. 3.3 and 3.4]. 



Appendix B. Global existence results for systems with time dependent / 



The main result of this section is Theorem |B.4| a special case of Theorem 3.8 It is applicable 



to solutions of Problem (2.8). Theorem B.4 is enough for our purposes as seen in the examples 
in SJ5l To prove the Theorem, we will modify the proof of Morgan 1989 Th. 2.2] with stronger 
control of the parameter r, that appears in (3.11). 



We start with two Lemmas from Morgan 



the results of Ladyzhenskaya et al. 1968 



1989| Lem. 4.1, Lem. 4.2, (4.12)] which follow from 



[Th. 9.1 p. 341]. We then use a duality approach to 



prove global existence of classical solutions to (2.8) 
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B.l. Lem ma (G lo bal ex iste nce) Let u be the solution of Problem (2.8). Let the f unct ion H fulfil 
conditions (3.8) — (3.10) in £3.6 and let the polynomial growth restriction on / (4.6) hold. Let 



Tm ax be as defined in (3.1) and suppose that, 



(B.l) 



for < t < T < T max and for all p e (1, . . . , oo) 
there exist M p , N p > and < 8 p < 1 such that 



E 



IIM-0|| Lp(OoX(T , T)) < M P (T - r) + N P (T r)\\H(n)\\l {nox{TT)y 



then T„ 



oo. 



B.2. Definition (Dual problem) A key ingredient of the proof of Theorem B.4 is the dual solution 
ip. Where for i = 1, . . . , to, ip is the solution of the scalar equation 



d t ip(€, t) = -DiAipfc t) - 6(£, t) for £ e fl and < t < T < T a 
( B.2 i { \ v ■ VV](€, *) = 0, I G dn Q , t e [0, T) 

^(S,T) = 0, |efi , 



where 6 > is such that, for all p € (1, . . . , oo), \\8\\l (n x[o,Tl) = 1- 



B.3. Lemma (Control of the solution to the dual problem (B.2 )) Let ip be as defined in B.2 For 
i = 1, . . . ,to and for pG (1, . . . , oo), there exists C p _t > such that, 



(B.3) 



l^ll Lp (no;L3o[0,T]) < C P,T 



We now state the main result of this section. Namely, the applicability of a special case of 



Theorem 3.8 to solutions of Problem (2. 



(2.8) admits a global classical solution, i.e., T lr 



oo (cf. (3.1)) 



B.4. Theorem (Sufficient conditions for global existence of solutions to Problem (2.8)) Let As- 
sumptions 2.3| 3.4 and 4.1 hold. Let H, f and r satisfy the conditions in ^3.6 with r < 1 (cf. 
(3.11)), i.e., Problem (3.4) admits a global classical solution by Theorem 3.8 Then, Problem 



Proof We proceed by contradiction. Assume T max < oo. Let u and / (cf. 4.2) be the solution 
and zero order ter m of Problem (2.8) respectively. From the proof of Lemma |4. 3 the polyn omia l 



growth restriction (4.6) is satisfied by / and H. Since T max < oo, Lemma 

k 



B.l 



implies that (B.l ) 



does not hold. Let j G [1, . . . , to] denote the smallest k for which 



i=l 



|/li(Wj)|| 



does 



L p (O x(r,T)) 

not satisfy (B.l I. From the proof of Lemma 4.3 the intermediate sum condition ( |4.5[ ) is satisfied 
for with r = 1 (cf. (4.5)). From Lemma A.l we have the a priori estimate (3.14). 



We will show (3.14) and (4.5) imply (B.l) is satisfied for j, obtaining a contradiction. We split 



the remainder of the proof into steps. 
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(B.4) 



(B.5) 



(B.6) 



(B.7) 



We first show the following inequality (corresponding to Morgan 1989 (4.6) — (4.9)]): 
For < T < T max 



r T 



Jn Jo l=1 u i 

Jo Jn a l=1 \ u j ) 

+ i>(£, s) ((fc x C? + k 7 )H (u(t s)) + fc 8 ) d s 

j 
+ $>^(t0)^(u°(O)d| 



--h+h + h- 



From ( B.2 ) we have for i £ [1, . . . , m] 



o Jn 

= 11 -h i (u i (Z,s))[d s i> + D j Ailj](Z,s)dSd. 
Jo Jn 

hi (ui(£,s))d s ip(£,s) 



< 



JQ 



D, 



xf}(^,s)D i A^h i (ui(^s))d^ds, 



where we have used integration by parts and the homogenous Neumann boundary con- 



ditions. From Problem (2.8| and the convexity of H (3.9|, we have 

-/ f ^,s)D i A^h i (u i ($,s))d$ds 

Jo Jn 



o Jn 

T 



o Jo, 



1>{£>*)K («< (t «)) /i(«(€, *),«) 



-/»<(«*(€, a)) fl^(€,«)d€da+ / ^(€,0)ft,- (u?®) d£, 

./fin V 7 



where we have used integration by parts and the final condition of (B.2). Combining 



(B.6) and (B.5), we obtain 



/ ( hi {ui{£, s))e&s) d$d s< f [ (^-lW^(£,s))^(£, 

Jo Jn Jct Jo K-Lsi J 



s) 



D, 



+ ^(Z,s)ti i (u i (Z,s))f i (u((-, S ), S )d S 



D, 



^(6 o)fc («?($)) <u. 
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Summing (B.7) over i < j and using condition (|4.5[) with f = 1, we obtain (B.4). For 



the case j = 1, we have introduced the convention /J(-) = 0. 



Step 2: We shall use Lemma B.3 and (B.4) to obtain the following inequality (as in Morgan 



1989 (4.10)— (4.16)]): For all p € (l,...,oo) there exists K VtT > and < 5 P < 1 



independent of u and 9 such that, 



(B.8) f f £a ji ^h i fa(t,8))efos)dSda<K p , T (l + \\H(u) 

Jo Jn i=1 u j \ 



L p (O x[0 : T]) 



Let p,q € (1, . . . , oo) be such that, - + - = 1. Dealing firstly with I\ (cf. (B.4)), wc 



have by Holder's inequality and the regularity estimate (B.3) 

D 



j-i 



'^E 



« ; 



(B.9) 



i-i 

2 = 1 



1 



D: 



Cq,TT\\hi (u % ) \\ L ^ noX[0 ^ 



1- 






c,, r r m p + jv p ||^(u) 



L p (O x[0 : T]) 



where we have the assumption that (B.l I is valid for i < j. Dealing with / 2j we have by 



Holder's inequality and the regularity estimate (|B.3) 
h< I HV&Ollwo.ri 



y ((A :i c 1 2 +fc 7 )ij(«(^, S ))+A :8 d S jdi 



(B.10) 



<C q , T (hCf + k 7 ) 



H(u(; S )) 



fcsT 



<C q , T (g{T)\Cl Q \ 1/p + k & T 



L q (Q ) 



for some g € C[0, oo). Where we have used the a priori estimate (3.14). Finally dealing 
with I 3 using Holder's inequality and estimate (B.3) we have 



(B.ll) 



J 

i=l 
j 

< 2_^ O'jiCq^Cp, 



L p (n x[0.T]) 



where we have used the boundedness of u° and condition (3.9). Combining (B.9), (B.10), 



(B.ll) and (B.4| yields (B.8). 



Step 3: We now show (B.l ) holds for j. Let p,q <E (1, . . . , oo), be such that, - + - = 1. We recall, 



from (3.9) and Definition B.2 that for i e [1, . . . , m], hi, 9 > and ||#||l (si x(o,t)) = 1- 
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Using duality we obtain 



<</ 



'3 i 






^II^(^)IIl p (o oX (0,T)) 

4 = 1 

. / dA 



(B.12) 



■fifr^D- 



f f Y / h i (u i ($,s))6(t,s)d$ds 



'o i=l 



< 



ft t* 






i(€,*))0(6*)d€d* 



<^,T(H-||F(ii)||- (fio)<[0jT]) 



where we have used (B.8|. 



Thus we have a contradiction and we conclude T max = oo completing the proof. 



□ 



The proof of Theorem 3.8 follows from more technical use of Holder's inequality in (B.lOl and 
in the case of assertion (3.17[) we also require the a priori estimate (3.15). We refer to Morgan 



1989 (4.13)— (4.16) and (4.27)— (4.19)] for specific details. 
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